rm(list = ls(all = TRUE))
library(foreign)

##############################
## Endogenous Primogeniture ##
##############################
d <- read.dta("KS_replication_european_monarchies.dta")
d$changed <- ifelse(d$country == "Bohemia" | d$country  == "Denmark" | d$country == "England" | d$country == "France" | d$country == "Gwynedd" | d$country == "Hungary"
                    | d$country == "United Provinces" | d$country == "Norway" | d$country == "Poland" | d$country == "Russia" | d$country == "Scotland" | d$country == "Sweden", 1, 0)
dsub <- d[d$primogeniture == 0 & d$defactoprimogeniture == 0, ]

## Deposal
pA <- mean(na.omit(dsub$dow_deposed[dsub$changed == 1]))
ndat <- dsub[is.na(dsub$dow_deposed) == FALSE, ]
nA <- nrow(ndat[ndat$changed == 1, ])
vA <- (pA * (1 - pA)) / nA
ci.upperA <- pA + qnorm(0.975) * sqrt(vA)
ci.lowerA <- pA - qnorm(0.975) * sqrt(vA)
datA <- c(ci.lowerA, pA, ci.upperA)

pB <- mean(na.omit(dsub$dow_deposed[dsub$changed == 0]))
nB <- nrow(ndat[ndat$changed == 0, ])
vB <- (pB * (1 - pB)) / nB
ci.upperB <- pB + qnorm(0.975) * sqrt(vB)
ci.lowerB <- pB - qnorm(0.975) * sqrt(vB)
datB <- c(ci.lowerB, pB, ci.upperB)

diffAB <- pB - pA
ci.upperAB <- diffAB + qnorm(0.975) * sqrt(vA + vB)
ci.lowerAB <- diffAB - qnorm(0.975) * sqrt(vA + vB)
datC <- c(ci.lowerAB, diffAB , ci.upperAB)

pdatB <- rbind(datA, datB, datC)

## Figure 1, left panel
pdf("figure1left.pdf")
par(mar = c(5, 5, 0, 0))
plot(c(1, 2, 3), pdatB[, 2], type = "n", xaxt = "n", xlab = "", ylab = "", xlim = c(0.7, 3.3), yaxt = "n", ylim = c(0, 0.6), frame = FALSE, cex.axis = 1.5)
points(c(1, 2, 3), pdatB[, 2], pch = 19, col = "black", cex = 1.5)
for(i in 1:3) lines(c(c(1, 2, 3)[i], c(1, 2, 3)[i]), c(pdatB[i, 1], pdatB[i, 3]), lwd = 1.5)
axis(1, at = c(1, 2, 3), c("Adopted\nPrimogeniture", "Never Adopted\nPrimogeniture", "Difference"), tick = FALSE, cex.axis = 1.5)
axis(2, at = seq(0, 0.6, 0.1), labels = seq(0, 0.6, 0.1), las = 2, cex.axis = 1.25)
mtext("Probability (and difference in probability) of Deposal", 2, cex = 1.5, outer = TRUE, at = 0.55, line = -1.25)
dev.off()
graphics.off()


## Figure 1, right panel
a <- density(dsub$tenure[dsub$change == 0])
b <- density(dsub$tenure[dsub$change == 1])
ks.test(dsub$tenure[dsub$change == 0], dsub$tenure[dsub$change == 1], alternative = "greater")
pdf("figure1right.pdf")
par(mar = c(5, 5, 0, 0))
plot(a, xlab = "", ylab = "", main = "", yaxt = "n", frame = FALSE, ylim = c(0, 0.1), cex.axis = 1.25)
lines(b)
polygon(a, col=rgb(0, 0, 1, 0.5), border=NA)
polygon(b, col=rgb(1, 0, 0, 0.5), border=NA)
text(36, 0.09, "D = 0.246\np-value = 0.000", pos = 4, cex = 1.5)
text(25, 0.06, "Never Adopted Primogeniture", cex = 1.5)
text(27, 0.03, "Adopted Primogeniture", cex = 1.5)
axis(2, at = seq(0, 0.10, by = 0.02), labels = seq(0, 0.10, by = 0.02), las = 2, cex.axis = 1.25)
mtext("Tenure (Years)", 1, cex = 1.5, at = 0.55, outer = TRUE, line = -1)
mtext("Density", 2, cex = 1.5, at = 0.55, outer = TRUE, line = -1.3)
dev.off()
graphics.off()
